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ABSTRACT 

The gravitational properties of a torus are investigated. It is shown that a torus can be 
formed from test particles orbiting in the gravitational field of a central mass. In this 
case, a toroidal distribution is achieved because of the significant spread of inclinations 
and eccentricities of the orbits. To investigate the self-gravity of the torus we consider 
the TV-body problem for a torus located in the gravitational field of a central mass. 
It is shown that in the equilibrium state the cross-section of the torus is oval with 
a Gaussian density distribution. The dependence of the obscuring efficiency on torus 
inclination is found. 
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1 INTRODUCTION 

' Obscuring dusty tori are important structural elements of 
[ Active Galactic Nuclei (AGN). In the framework of the 
. unified scheme, the difi^erences in the AGN are explained 
' by the orientation of the torus relative to an observer 
, (Antonucci 1993; Urry & Padovani 1995). Full obscuration 
■ of the central engine (black hole + accretion disk) and the 
' Broad Line Region (BLR) is realized when the torus is 
, seen edge-on. In this case, only narrow emission lines are 
> detected in spectrum, and this type of AGN is called 'type 
' 2'. In contrast, if the torus is inclined at some angle to the 
, line of sight, the emission from the accretion disk and BLR 
can be seen to the observer ('type 1'). Following from the 
statistical data, the torus must be geometrically thick in 
order to explain the observed properties of Seyfert galaxies 
(Osterbrock & Shaw 1988; Schmitt et al. 2001). 

The first direct evidence for the existence of a compact 
structure in the nucleus of Sy2 galaxy NGC 1068 was found 
with help of near-infrared bispectrum speckle interferometry 
(Wittkowski et al. 1998; Weigefi et al. 2004) and speckle 
interferometry (Weinberger, Neugebauer & Matthews 1999). 
Subsequently, observations of NGC 1068 using VLTI/MIDI 
in the infrared band (lO/^m) allowed, for the first time, a 
spatial distribution of temperature in the torus to be obtain 
(Jaffe et al. 2004). It was discovered that the torus has two 
components: a hot compact component (T > 800A') and 
a warm component (T = 300 A). The hot component is 
about 1.35 parsec long and 0.45 parsec thick; the size of the 
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warm component is 3 x 4pc (Raban et al. 2009). Thus, the 
observations confirmed a substantial geometrical thickness 
of the torus: the minor-to-major radius ratio ro is about 
0.7. It is suggested that the hot component discovered by 
JafFe et al. (2004) is the inner funnel of the obscuring torus 
heated by radiation of the accretion disc, while the warm 
component is the torus body itself (Schartmann et al. 2005, 
DuUemond & Bemmel 2005, Raban et al. 2009). 

The mass of the obscuring torus can be estimated from 
the rotation curve obtained from observations of maser 
emission. Greenhill et al. (1996) obtained VLBI synthesis 
images of the H2O maser emission from the central region 
of NGC 1068. They found a good fit to the data by assuming 
a non-Keplerian rotation curve of the form V oc r""'^^ on 
spatial scales above 0.6 pc. The deviation of the rotation 
curve from Keplerian may be a result of the infiuence of the 
self-gravitation of the torus (Greenhill et al. 1996). It turns 
out that in the case of NGC 1068, the torus mass could be 
comparable to the mass of the central black hole (Hure 2002; 
Lodato & Bertin 2003). 

One of the problems concerning obscuring tori is how 
they stand up against gravity. To explain the required 
geometrical thickness of a torus, the magnitude of random 
speeds in the vertical direction must be about: AI4/Vorb ~ 
ro. Krolik & Begelman (1988) noted that if this motion is 
thermal, the corresponding temperature (~ lO^A) would 
be far too hot for dust to survive. Therefore it was 
assumed that the torus material is distributed in clouds. 
To support cloud motions against losses by collisions, 
orbital shear energy can be randomized if magnetic fields 
make the clouds sufficiently elastic (Krolik & Begelman 
1988). However, the required magnetic fields have not been 
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detected. Pier & Krolik (1992) suggested a mechanism in 
which a geometrically thick torus can be explained by the 
influence of the infrared radiation pressure. Krolik (2007) 
obtained an explicit solution of hydrostatic equilibrium 
within a torus on the assumption of a continuous medium. 
In the framework of this solution it was shown that the 
infrared radiation pressure can support geometrically thick 
structures in AGN. In this case, the density distribution 
in torus cross-section reaches its maximum near the inner 
edge of the torus with a subsequent decrease to the outer 
edge. Schartmann et al. (2005, 2008) suggested a turbulent 
torus model in which the dusty clouds are ejected by stars 
and take over the stellar velocity dispersion. In this model, 
the torus is considered as a continuous medium and the 
turbulent motions lead only to a weak stabilization against 
gravity (Schartmann et al. 2010). Starburst-driven tori were 
investigated with the help of a hydrodynamical simulation 
of the interstellar medium by Wada & Norman (2002) and 
Wada, Papadopoulos & Spaans (2009). The other idea is 
the presence of a dipole vortex motion in the torus that can 
support its geometrical thickness (Bannikova & Kontorovich 
2007). There are also models in which the obscuring region 
is produced by a hydromagnetic wind (Konigl & Kartje 
1994; Elitzur & Shlosman 2006). The main problem of these 
models is the unknown origin of the large-scale magnetic 
field needed to support such winds. An alternative to models 
of wind scenarios is that of outflows driven by infrared 
radiation pressure on dust (Dorodnitsyn, Bisnovatyi-Kogan 
& Kallman 2011, 2012). 

In hydrodynamic models the torus is considered as a 
continuous medium, and hence the motion of matter (the 
clouds) is parallel to the torus plane of symmetry. In this 
case, gravity tends to compress the torus into a disk, and 
therefore the substantial forces are required to work against 
gravity. In this paper we consider discrete clouds moving 
in the gravitational field of the central mass. In this case, 
the orbital plane of each cloud passes through the central 
mass. The clouds form a Keplerian disc if their orbits lie in 
one plane. The toroidal structure can be formed by clouds 
moving in orbits of different inclinations and eccentricities. 
Such orbits are similar to those of stars moving around the 
supermassive black hole in the Galactic centre (Ghez et al. 
2005). Note that recent observations have detected in this 
region a dusty cloud, which also moves in a very elongated 
orbit with large inclination (Gillessen et al. 2012). 

In this paper we investigate the gravitational properties 
of a torus composed of clouds moving in orbits with different 
inclinations and eccentricities. We start with an idealized 
case, considering clouds in approximation of test particles 
moving on Keplerian orbits in the gravitational field of the 
central mass (Section 2). In Section 3 we considered the 
problem of test particle motion in the inner gravitational 
potential of a homogeneous circular torus and central mass. 
In Section 4, the A'^-body problem for a torus located in 
the gravitational field of a central mass is investigated. The 
equilibrium cross-section of self-gravitating torus is found, 
and its evolution over long times is analyzed. The obtained 
results are used for the interpretation of the observed 
features of the obscuring tori in AGN (Section 5). 
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Figure 1. The scheme of the torus, which consists of particles 
orbiting in gravitational field of a central mass Mc- 

2 KEPLERIAN TORUS 

Let us consider the problem of test particles motion in the 
gravitational field of a central mass A'lc (Fig.l). This is a two- 
body problem for each particle. We can impose conditions on 
the orbital elements such that the particles form a toroidal 
structure. In this case, a substantial geometrical thickness of 
the torus can be reached owing to the significant scatter of 
inclinations and eccentricities of the particles orbits. Because 
particles move on unperturbed (Keplerian) orbits, we will 
call a torus formed in this way a "Keplerian torus". 

To create a Keplerian torus we impose the following 
conditions: 

a) the shape of the cross-section of the torus must be 
close to circular or elliptic; 

b) the spatial distribution of particles in the torus must 
be close to homogeneous; 

c) the orbits of particles in the comoving system of 
coordinates must be nested. This condition allows us to 
avoid intersections of orbits and hence collisions between 
particlesQ 

The torus is characterized by two geometrical 
parameters: major R and minor _Ro radii (Fig. 1). It is 
convenient to introduce the geometrical parameter ro = 
Ro/R. To ensure a quasi-circular cross-section of the torus, 
we should choose the semi-major axes of the orbits of 
all particles equal to the major radius of the torus. In 
order to keep all the particles under the torus surface, the 
eccentricities of their orbits must be in the range e — [0, ro] 
and the inclinations, measured from the symmetry plane of 
the torus, must satisfy the following condition: 

i = arcsin | q — , | . (1) 

where g is a parameter allowing us to change the ellipticity 
of the torus cross-section. Let p = ^/x^ + y^ / R be the radial 
coordinate of a particle in the torus plane of symmetry, 
measured from the central mass and normalized to R. Then 
the coordinates of the inner and outer boundaries of the 

^ Note that in Section 4 we will pass from test particles to more 
realistic objects that have sizes and masses. Therefore collisions 
between particles will be taken into account. 



torus are determined by the eccentricity of the orbit pmm = 
1 — e and pmax = 1 + e. Let us suppose that the argument 
of perihelion for all orbits oj = 0; that is, the line of nodes 
coincides with the semi-major axis. We generate randomly 
the longitude of the ascending node and the true anomaly 
V in the range from to 2ti. Knowing the orbit elements, we 
can calculate spatial coordinates for each particle using the 
well-known formulae (Duboshin 1968) 



(x,?/, 2) = r X (a,/3,7), (2) 



where 



a = cos Q. cos V — sin Q, sin v cos i 
P = sin Q cos 1/ + cos Q sin u cos i 
7 = sin ly sin i 

and the module of the radius vector from the central mass 
to the particle 



(3) 



1 + e cos u 

The velocity components of the particles are determined in 
the following way: 





(1 + e cos 



(4) 
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e ) is a focal parameter, and / — y^GMc ■ p is a 
module of the kinetic moment. Thereafter, we will use the 
system in which G — R = Mc = 1. Since R — 1, the value 
of the geometrical parameter ro uniquely determines the 
minor radius of the torus and consequently its geometrical 
thickness. Finally, an algorithm to simulate Keplerian torus 
is the following. 

1. Set the number of particles A'^, the geometrical 
parameter ro and the ellipticity of the torus cross-section 

2. Generate randomly the eccentricities within the given 
limits e — [0, ro] for each particle and find corresponding 
values for the inclinations by the formula ([T]). Also generate 

and u randomly for these orbits. 

3. Find the spatial coordinates and velocity components 
by elements of the orbits according to - (SI) and use them 
as the initial values. 

4. Solve the equations of motion with the obtained 
initial conditions 0. 

Because the investigated system is symmetric with 
respect to azimuthal angle, it is convenient to use a comoving 
system of coordinates. This system of coordinates is formed 
by a plane perpendicular to the plane of torus symmetry 
and passing through the radius vector of a particle. Thus, 
one axis of this system p coincides with a projection of the 
radius vector of a particle on the plane of torus symmetry, 
while second one ( — z/R coincides with the axis of torus 
symmetry. Because a particle is orbiting around the central 



^ The results in the form of animations are presented at 
|http: / /astrodata. univer.kharkov.ua/ agn/torus/| 
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Figure 2. Trajectories of particles with the eccentricities e = O.lfc 
(fc = 1..7) in a comoving system of coordinates: a) g = 1, b) 
5 = 0.7. The central mass is located at the point with coordinates 
(0,0). Dashed lines show the circular (a) and elliptical (b) cross- 
sections. For the limiting case of a thick torus (outer solid curves), 
the cross-sections of a Keplerian torus differ noticeably from the 
corresponding circular (a) or elliptical (b) ones. 



mass, the comoving system is also orbiting with a velocity 
equal to the orbital velocity of the particle. The comoving 
system of coordinates is convenient because the trajectories 
of particles in this system reflect the shape of the torus cross- 
section. 

The equation of a particle trajectory in the comoving 
system of coordinates for the case of a Keplerian torus has 
the following parametric form: 



P = 



1 + e cos p 



1 - g2 



1 - 



( = qe\/l — e2 



1 -I- e cos u 



(5) 



(6) 



The trajectories of the particles in the comoving system of 
coordinates calculated with ((Sj - ((61 are presented in Fig 2. 
The shape of each orbit determines the shape of the cross- 
section of the torus with a given value of ro. The period 
of a particle in comoving system of coordinates equals the 
orbital period. Furthermore, the periods of all particles are 
equal because the semi-major axes of all orbits equal unity 
{R=l). 

It is seen from Fig.2 that the cross-section of the 
Keplerian torus is different from the circular cross-section 
for large values of orbit inclinations (the outer trajectories 
in Fig.2). At small inclinations, the trajectories are 
nearly circular (inner curves in Fig. 2a). Interestingly, the 
trajectories at large inclinations make the cross-section 
of the torus a more complicated shape with two humps 
(Fig.2a). The restriction on the geometrical thickness of a 
Keplerian torus follows from Indeed, because sini < 1 
we obtain the ma:ximum value of the eccentricity emax ~ 
l/yg^ + l from {!]) which is the upper restriction on ro. 
For q = 1, the highest possible value of the geometrical 
parameter of a Keplerian torus is ro — l/\/2 « 0.7. This 
limit corresponds to the maximum value of orbit inclination 
i = n/2; that is, the outer particle is moving in the plane 
perpendicular to the torus plane of symmetry in this case. 
Note that if the eccentricities and the inclinations of the 
orbits of all particles tend to zero, then the Keplerian torus 
degenerates to an infinitely thin ring. Obviously, because 
this torus is formed from test particles (interaction between 
particles is not accounted for), the orbits of the particles are 
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always closed and the Keplerian torus is stablelf] This is a 
limiting case for a self-gravitating torus when its mass tends 
to zero. In reality, the torus has a mass, and, consequently, 
one needs to account for its gravitational potential. In the 
next section we will investigate the trajectories of a test 
particle in the gravitational potential of the central mass 
and a homogeneous circular torus. 



3 TRAJECTORIES OF TEST PARTICLE IN 
THE INNER POTENTIAL OF A 
HOMOGENEOUS CIRCULAR TORUS AND 
THE CENTRAL MASS 

Consider the motion of a test particle in the inner 
gravitational potential of both the torus and the central 
mass. Because the expression for the torus potential is not 
integrable explicitely, we can take advantage of its expansion 
for the analysis of particle trajectories. Bannikova, Vakulik 
& Shulga (2011) obtained an approximate expression for the 
inner gravitational potential of the homogeneous circular 
torus in the form of a power series up to second-order terms: 



inner / > \ 
Vtorus (??,<,; ''Oj 



GMtorv 

2tvR 



c + ai^ + a2( — ] +b2( — 



(7) 

where Mtoms is the torus mass, 77, C, are dimensionless 
coordinates [j] = p — 1) and the coefficients depend on the 
geometrical parameter tq: 



ai = 8A:(1 + lnfc), k = ro/8 
a2 = -1 - 4fc2(ll + lOlnfc), 62 



-l + Ak\3 + 2lnk). 



The approximate expression ^ describes the inner potential 
of a homogeneous circular torus with good accuracy up to 
ro = 0.5 and allows to perform a qualitative analysis of 
the trajectories in the gravitational system "torus + central 
mass". The coefficient ai defines a shift of the potential 
maximum with respect to the centre of the torus cross- 
section and the coefficients 02, b2 are related to the deviation 
of the torus potential from quadratic law. As noted in a 
previous work (Bannikova, Vakulik & Shulga 2011), the 
inner potential of the torus can be represented as the sum of 
a cylinder potential and a term comprising the curvature of 
the torus surface ip{r) = (pcyi + fcurv Two types of central 
fields are known in which all trajectories of finite motions 
are closed: (p{r) oc 1/r and (p (x (Landau & Lifshitz 
1988). Obviously, closed orbits are possible in the cylinder 
because its gravitational potential is ipcyi oc r^. However, the 
presence of a potential of curvature tpcurv brings to forming 
open and more complex orbits of the particle in the inner 
potential of the torus. We investigate this problem in detail. 
For the qualitative analysis of particle trajectories in the 
inner region of the torus, we restrict ourselves to the case of 
a torus with circular cross-section and homogeneous density 
distribution. Since the potential of the considered system 
is axisymmetric, the projection of the angular moment of 
the particle on axis is a constant value (/,; = Const). In 



^ Because there is an analitical solution for r(i') in case of a 
Keplerian torus we use it as the initial condition in the A'^-body 



this case, it is convenient to introduce an effective potential 
(Binney & Tremaine 1994): 
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where the potential of the central mass Mc 
GM,, 1 



(8) 



(9) 



We can represent v^e// in the form of a power series up to 
second-order terms: 
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where the coefficients 

£ = 27r/i(2 - Z^) + c 

di = 2iTfiro{2 — Y) + ai 

d2 = — 7r/ir-o(3/^ — 2) -f- a2 

0,2 = —TVfirQ + &2 

= Mc/ Mtorus is the ratio of the central mass to the 
torus mass, and the dimensionless parameter I defines a 
fraction of the angular momentum from the Keplerian value 
— YGMc/R. It was shown in (Bannikova, Vakulik 
& Shulga 2011) that the maximum of the potential of 
the homogeneous circular torus is displaced inside from 
the centre of its cross-section. As a result, the torus 
must be compressed along the major radius. To prevent 
this compression, an orbital motion is necessary: the 
gravitational force tending to compress the torus along the 
major radius is compensated by the centrifugal force. This 
corresponds to a shift of the maximum of the effective 
potential to the centre of the torus cross-section; that is, 
fii = in (|10p . whence it follows a condition on the 
coefficient of the angular momentum: 



I = 1 — ai/(27r/iro). 



(11) 



The equations of a particle trajectory in the comoving 
system of coordinates can be obtained by solving the 
equations of motion: 

d^C _ 1 d(peff 



d rj _ 1 difieff , 



dt^ R dr) ' ^ dt^ R d( ' 
These equations have the form of a harmonic oscillator 



(12) 



(13) 

Then the solution of the equation (|13p can be presented in 
parametric form: 



■q{t) = rio cosftJA ■ t + a], 
C{t) = (ocoslujB ■t + P], 



(14) 



where 770, Co are the coordinates of the particle at initial 
time, and 

2 GMtorua r 2/r,72 o^ 1 /1t;^ 

= _ nc,^2 [t^''o(3/ - 2) - a2\ , (15) 



nR^r 



2 GMtor 



\nfiro 



(16) 



simulation (Section 4). 



Thus, the coordinates of the particle (|14|) evolve like the 
displacements of two harmonic oscillators with the epicycle 
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Figure 3. Trajectories of a test particle in tlie inner potential of 
liomogeneous circular torus and the potential of the central mass. 
Left column: the trajectory in the comoving system of coordinate 
for a) = 0.4; I = 0.95, c) = 0.7; I = 0.90. The centre 
of the torus cross-section is marked with a point. Right column: 
projection of the trajectory on the torus plane of symmetry for b) 
V( = 0.4; I = 0.95; d) = 0.7; I = 0.90. Dashed Unes correspond 
to edges of the torus cross-section. For all cases fi = I, ro = 0.5, 
G = R = 1, t = 150. 



to A and vertical ojb frequencies, respectively. In this case 
the orbit is open and the trajectory in the comoving system 
of coordinates fills some rectangular region during a period 
Pbox whereupon this process is repeated. These orbits are 
known as box orbits (Binney &Tremaine 1994). It is seen 
from (jlip that for firo ^ 1 the coefficient of angular 
momentum I — >■ 1 and the ratio of epicycle frequency uja 
to vertical frequency ljb takes the form: 



UJA _ 




- a2 


UJB 


1/ ■K^n■l 


-b2 



(17) 



Then, the period of a box orbit, expressed in orbital periods, 
is determined by the expression Pbox = 1/|1 — f\- For ro — 
0.3 and fj. = 50; 25; 1/0.06 we obtain Ptox = 213; 113; 80. In 
the case of /iro ^ 1 (Z — >■ 1) it is apparent from p7p that 
iUA UJB- This limiting case corresponds to the case of a 
Keplerian torus. 

Fig. 3 presents the trajectories of a test particle in the 
inner gravitational potential of the homogeneous circular 
torus and the potential of the central mass. For this we used 
an expansion for the inner potential of the torus in a power 
series up to fourth-order terms: 



znner / /- \ 
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(Bannikova, Vakulik & Shulga 2011). In contrast to a square- 
law potential ((7|, in this case we consider a more explicit 
approximation of the inner potential of the torus fTS]). 
Therefore, the form of the region bounding the place of the 
particle in the torus cross-section has more complex shape 
than (|14p (Fig.3a,c). Note that box orbits arise as a result 
of mistiming of the orbital period of the particle and its 
period in the comoving system of coordinates. If the torus 
mass is increased, the orbital period of particles is decreased 
due to turning the line of nodes. This leads to a shift of the 
pericentre and, therefore, to the formation of a rosette orbit 
in the torus plane of symmetry (Fig.3b,d). When the torus 
mass is decreased, the trajectory in the comoving system of 
the coordinates fills the box area for a longer time; that is, 
the period of the box orbit increased in this case. 



4 SELF-GRAVITATING TORUS IN THE 

GRAVITATIONAL FIELD OF A CENTRAL 
MASS: iV-BODY SIMULATION 

This section is devoted to an investigation of a self- 
gravitating torus located in the field of a central mass. We 
restrict ourselves to the case of a torus whose mass is up 
to 10 per cent of the central mass. As the initial condition 
we use a Keplerian torus (Section 2). We investigate the 
following main questions: 

1. Are the tori stable on time scales comparable to the 
lifetimes of astrophysical objects? 

2. What is the shape of the cross-section of a torus in 
an equilibrium state? 

3. What is the density distribution of particles in the 
equilibrium state of the torus? 

As will be seen below, the self-gravitation significantly 
changes the cross-section shape of the torus as compared 
with the case of a Keplerian torus. The gravitational 
potential of a torus with mass Mtorus, composed of N 
gravitating particles of equal masses, at an arbitrary point 
r(p. A, (") has the form 



G Alt or 

w 



■E- 



1 



(19) 



where r{p, X, is the dimensionless radius vector of i-th 
particle normalized to the major radius of the torus R. 
Taking into account the gravitational interaction of particles 
forming the torus, the forces acting on a particle can be 
represented as the sum of regular and irregular components. 
The regular component is associated with a 'smoothed' 
potential, and irregular one appears because of the direct 
gravitational interaction between nearest particles. The role 
of the irregular forces increases in the case of a relatively 
small number of the particles. Therefore, the torus potential 
can be represented as follows: 



'Ptorus{p,X,0 = ftoLsipX) + 'PTorus{P-A,C), 



(20) 



where the coefficients of the expansion can be found in 



where the first term is the regular part of the torus potential 
that is independent on the azimuth angle A, while the second 
term is an irregular (random) part of the potential defined 
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by 



R 



^^torus / 

^-^ r — 1 



K(r')dV' 



(21) 

Here K,{r) is the volume density of the torus. Accounting for 
the central mass leads to the following expression for the full 
potential: 



'^'{P, K C) = 'fic{P, C) + ftorusip, A, C), 



(22) 



where the potential of the central mass is defined by 
expression The limiting case ^plorua ~^ (neglecting 
interactions between particles) corresponds to the problem 
of the motion of a test particle in the regular gravitational 
potential of the torus. This problem was analyzed in detail in 
the approximation of a homogeneous circular torus (Section 
3), where it was shown that the trajectories of particles in 
the cross-section of the torus relates to the type of box 
orbits. In the following, we will show that the role of the 
regular part of the potential appears in changing the sizes of 
the torus cross-section for A'^-body simulation. On average 
(by averaging over all i-indices), ipYorus ~^ 0, but at any 
particular point iplorua can vary considerably from zero, 
taking both positive and negative values. As a consequence, 
the velocity of particles moving in the regular potential of 
the torus will be subject to random perturbations owing to 
the forces determined by the irregular part of the potential. 
As a measure of the irregular part of the potential we 
choose a potential created by a particle at a distance I — 
\/Vtorua/N equal to the average distance between particles: 
ipi — Grrii/l. The masses of all particles are the same 
{rrii = Mtorus/N), and therefore ipi oc N~^^^; that is, 
the irregular forces become significant for a relatively small 
number of particles N. 

The considered A''-body problem is reduced to the 
numerical integration of the equations of motion taking 
account of the central mass: 

GMc ri Fi 
-^-—^^ + - (23) 

where ai is the acceleration of i-th particle. The total 
gravitational force acting on i-th particle is 



Grrii 



ri 



(24) 



where £ is a parameter allowing us to avoid the infinite 
increase of the force of interaction between two particles 
when they come close together (Aarseth 1963, 2003). 
However, the introduction of this parameter means that 
we simulate not point particles but the mega-particles of 
spherical form with dimensionless radius e in which the 
density distribution obeys Plummer's law (Plummer 1911; 
Saslaw 1985). In the following numerical simulations we take 
the parameter e = 0.01. Note that the clouds penetrate 
through each other for collisions in this case. 

A considerable amount of computer time is needed for 
the numerical simulation of a system consisting of large 
numbers of particles and accounting for their gravitational 
interaction. In this research we used the technology of 
parallel computing CUDA developed by Nvidia, which 
allows us to significantly reduce the computation time by 
using the computing resources of graphics processing units 



(CPUs). This technology is especially effective for A^-body 
problems (Belleman, Bedorf & Portegies Zwart 2008). The 
numerical simulation was performed using the graphics 
card NVIDIA GeForce GTS 250 (192 CUDA Cores). The 
calculation time was about 2 hours for a system of 8 192 
particles orbiting over 1000 periods with time step equal 
1/150 of the period. 

Our main idea is that a torus should be more stable 
if it is formed by clouds moving in inclined orbits around 
the central mass. Therefore we use the Keplerian torus 
of the given geometric parameter ro (Section 2) as the 
initial condition. The solution for a Keplerian torus allows 
us to calculate the initial coordinates and velocities of 
particles with inclined orbits. The absence of quasi-periodic 
fluctuations of macro-parameters of a statistical system is 
a necessary condition for a system of particles to achieve 
a quasi-steady state. In order to test this condition, we 
control the values of the kinetic energy, potential energy and 
total energy of the system a few times over the period. In 
addition, in order to control the dynamics of the system, the 
statistical macro-parameters of the system of particles were 
recorded with the same frequency, namely the coordinates 
of the centre of the torus mass in its cross-section {rjc, Cc) as 
well as the effective sizes of the torus cross-section, which 
were calculated by the following expressions 



1 



1 



(25) 

where r^, are respectively the horizontal and the vertical 
sizes and the 'average' effective size of the cross-section: 



(26) 



Note that in the case of a homogeneous circular disc all 
the effective sizes, defined in such a way, are equal to the 
geometrical radius of the disc. The angular momentum of the 
system of particles in the comoving system of coordinates is 
defined by the expression 



(27) 



where the dot denotes a time derivative. Hereafter, this 
angular momentum (|27|) is referred as the 'fictitious' 
angular momentum. An analysis of the behavior of the 
fictitious momentum allows us to determine the degree 
of order (or chaos) of directions of particles movement 
in the torus cross-section. Furthermore, we periodically 
register the position (coordinates) of all particles in the 
plane of the torus cross-section. The accumulation of 
these positions during each 10 orbital periods allowed 
us to obtain the average density distribution of particles 
in the torus cross-section and to analyze the change of 
cross-section shape during the simulation. The evolution 
of the torus cross-section during the simulation for 
different ro is presented in the form of an animation at 
http://astrodata.univer.kharkov.ua/agn/torus' 



The simulation shows that large changes of kinetic and 
potential energy take place at the initial stage, t < 2Q 
(Fig.4), with the total energy conserved. The change of 
potential energy is related to the change of shape of the 
torus cross-section and to the density of particles in it. Just 
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Figure 4. The sum of twice the kinetic and the potential energy 
of system, 2Ei,i„ + Epot, during the first 100 orbital periods for 
Mtorus/Mc = 0.02 and N = 8 192. Initial conditions correspond 
to a Keplerian torus with ro = 0.3. 
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Figure 5. Variation of cross-section sizes of a torus with mass 
Mtorus = 0.02Mc, consisting of = 8 192 particles, during 1 800 
orbital periods. The initial conditions correspond to a Keplerian 
torus with ro = 0.3. 



in this short initial time interval, the distribution of particles 
in the torus varies significantly, and the torus acquires an 
equilibrium cross-section of oval shape. Furthermore, the 
sharper part of the oval cross-section is directed to the 
central mass, and the more flattened part to the outer region. 
This is the opposite to the situation in the initial state 
(Keplerian torus). Obviously, the significant fiuctuations in 
energy occur because of the tendency of the torus to acquire 
a cross-section shape inverted with respect to the initial one. 
For t > 20, the fluctuations of the kinetic and potential 
energy are negligible. This means that the shape of the torus 
cross-section does not change significantly anymore. We 
define this state in which the shape of the torus cross-section 
is no longer changing as the equilibrium state. Fig. 4 shows 
that the virial theorem is satisfied, namely 2Ekin + Epot ~ 0, 
which characterizes the stationary or linearly non-stationary 
state of the system. Fig.5 shows that the torus cross-section 
gradually increases with time. The gradual spread of the 
cross-section is related to the irregular forces arising from 
the gravitational interaction between particles. Indeed, the 
irregular forces dominate over the regular ones caused by 
the smoothed potential for the considered case of a 'light' 
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Figure 6. Variation of the fictitious angular momentum for a 
torus with mass Mtorua = 0.02Mc and N = 8 192 during 
1 800 orbital periods. Ig is normalized to the initial value that 
corresponds to a Keplerian torus with ro = 0.3. 



torus [Mtorus/Mc = 0.02). The periodic fluctuations in size 
(Fig.5) are related to the regular potential of the torus. It 
was shown in Section 3 that the orbit of a test particle 
in the inner potential of a circular torus in the comoving 
coordinate system is of the box-orbit type. In the case of A^- 
particles the box orbits are synchronized at the initial stage, 
which causes periodic fluctuations of the cross-section of the 
torus (Fig.5). The amplitude of these oscillations decreases 
at later stages because of the mistiming of the box orbits. 
As a consequence, the flctitious angular momentum in the 
comoving coordinate system, Iq, also decreases (Fig. 6). The 
system tends to a state with Ig — >■ 0, that is to the 
randomization the movement directions of particles in the 
torus cross-section. Note that the total angular moment of 
the system is conserved. The changes of sizes of the torus 
cross-section for a larger value of the initial geometrical 
parameter, ro = 0.7, are presented in Fig. 7. The significant 
difference from the previous case of a torus with the smaller 
initial value ro = 0.3 is obvious (Fig.5). Namely, the sizes of 
the cross-section change considerably during the initial time, 
which is connected to a significant change in the shape of 
the torus cross-section. Subsequently, the amplitude of the 
fiuctuations decreases and finally disappears after t = 300. 
This behavior can be explained as follows. In the case of a 
thick torus, the periods of box orbits in the initial stage are 
very different, which leads to a rapid mistiming of orbits in 
the torus cross-section. This confirms the analysis of change 
in the fictitious angular momentum (Fig.8). Thus, the cross- 
section of a self-gravitating torus spreads slowly (here with 
the vertical size increasing more slowly than the horizontal 
one). In the equilibrium state, the directions of particle 
movements in the torus cross-section are randomized. Note 
that if the central mass Mc — lO^M© then one orbital 
period [t — 1) in the centre of the cross-section of a torus 
with the major radius R — 2pc correspond to an average 
interval of 80 000 years (see Section 5). Thus, the results of 
numerical simulation for N = 8 192 (Figs. 5-8) for a time 
interval of 1 700 orbital periods correspond to 136 million 
years. 

When particles interact with each other, some of them 
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Figure 7. Variation of cross-section sizes of a torus with mass 
Mtorus = 0.109Mc, consisting of Af = 8 192 particles, during 
1 800 orbital periods. The initial conditions correspond to a 
Keplerian torus with ro = 0.7. 
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Figure 8. Variation of fictitious angular moment for a torus 
with mass Mtorus = 0.109Mc and N = 8192 during 1800 orbital 
periods. Ig is normalized to the initial value that corresponds to 
a Keplerian torus with ro = 0.7. 



lose energy and move to lower orbits, while the energy 
of other particles increases and they shift to more distant 
orbits or even escape the system. As a result of such 
interactions, the outer boundary of the torus cross-section 
gradually increases. By virtue of the random nature of 
particle interactions we can assume that the modulus of the 
total change of a particle velocity AV" is related to velocity 
changes in separate encounters of particles AVi as 



{AVf = ^(AFO 



(28) 



In the theory of stellar dynamics there is well-known 
expression for an estimate of the velocity dispersion {AVy 
reached during the interval At (Kulikovskii 1978, Saslaw 
1985): 



(29) 



(AV)^ = 32nGm^At • D ( — ) In 



B 



where m is mass of a particle, D is number of particles 
per unit volume, which for a homogeneous distribution is 
proportional to A^. Omitting the constant coefficients and 
assuming that the mass of particles m = Mtoms/N , the 
expression (|29|l takes the form 



(AVy oc G- 



N 



-At In 



B 



(30) 



As an estimate of the maximum impact parameter hmax, 
the average distance I between particles is commonly used. 
Then, the expression for the impact parameter is 



27r2i?3rg 



N 



and 



B 



(31) 



(32) 



In the considered dimensionless system we obtain the 
estimate 



In 



In 



iV2/3 



(33) 



Neglecting the weak logarithmic dependence (AV)'^ of A'' 
and Mtorus we finally obtain 



(34) 



Although the expression (|34|l was obtained for a simply 
connected system of particles, it works satisfactorily in the 
case of a torus, which is confirmed by results of simulation. 
For example, if A*'i is number of particles in the torus 
and Ati is interval during which a certain value of the 
velocity dispersion is reached, then it follows from (|34p 
that At2 — Ati ■ N2/N1; that is, the same value of the 
velocity dispersion for N2 = N\/2 will be reached during 
Ati = Afi/2. Variations of cross-section sizes for tori 
with different numbers of particles are shown in Fig. 9. 
Black curves were obtained from simulation of a torus with 
iVi = 16 384 for the interval ti = 3000. The gray curves 
correspond to the torus that consist of N2 — Ni/2 = 8 192 
particles for interval At2 = Ati/2 = 1 500 with subsequent 
expansion of the obtained curves by a factor of two. The 
black and gray curves coincide well, which indicates that 
the expression (|34p is satisfactory. Thus, provided we obtain 
results of the simulation of the system with a smaller number 
of particles in a relatively short time interval, expression 
^34\ ) allows us to predict the state of the system consisting 
of a larger number of particles for larger time intervals. In 
a similar way, the condition of re-scaling \34^ can be used 
to investigate systems with different values of Mtorus ■ 

We discovered from simulations that in the equilibrium 
state the torus has an oval cross-section with a density 
distribution that changes exponentially. The stability of the 
torus is achieved because the clouds move on inclined orbits. 
The distribution of particles in the torus cross-section was 
obtained by averaging their coordinates on the interval of 
100 orbital periods with a subsequent approximation by a 
function of the form 



n{vX;ro) = no(ro)exp[-/(7),C;ro)], 



(35) 



where no is the number of particles at a maximum of 
the distribution. For simplicity we set no = 1- Obviously, 
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Figure 9. Variation of cross-section sizes of a torus with initial 
conditions ro = 0.3, Mtorus = 0.02A'Ic- The black curves are 
obtained by simulation with A^i = 16 384 for the interval A4i = 
3000. The gray curves correspond to case N2 = Ni/2 = 8 192 at 
At2 = 1 500 with subsequent expansion of the obtained curves 
by factor of two. 



A steady state of the clouds can exist if the rates of the 
two processes are comparable (Krolik & Begelman 1988). 
Our point is that on average these processes do not change 
considerably the spread of orbits inclinations. Wind could 
be an additional source of clouds moving on such orbits: the 
clumpy wind coming off the accretion disk can feed the torus 
(Nenkova et al. 2008 Part II). In this case such clouds could 
acquire inclined orbits. 



5 APPLICATION TO THE OBSCURING TORI 
IN AGN 

The obtained results (cross-section shape and density 
distribution) can be directly applied to the obscuring tori 
of AGN. We will use the condition of obscuration obtained 
from analysis of the spectral energy distribution: for total 
obscuration of the central engine, the number of clouds in 
the equatorial plane must be greater than 5 (Nenkova et al. 
2008 Part I,II). 



ln(n/no) = / (»7, C) is equation of an isodence at level n/no 
for a given density distribution (|35p . Because the isodences 
have an oval shape we approximate the function / as a power 
series in the coordinates: 



(36) 



Fig. 10 shows an average density distributions in the 
equilibrium cross-section of a torus consisting of A'^ = 8 192 
particles for various values of the torus mass and various 
initial ro. These cross-sections were achieved in 1000 orbital 
periods. The parameters {Mtorus and ro) in Fig. 10 are 
chosen so that in all cases the tori have the same values 
of the initial volume density: utorue = Mtorus / (^tt^ R^Tq) = 
Const. We also performed special simulations by changing 
the initial shape of the torus cross-section. The results of 
the simulations showed that changes of the initial shape of 
the torus cross-section (the ellipticity parameter g in ([T])) 
does not substantially influence the equilibrium shape of the 
cross- sect ion. 

So, we come to the following conclusion. 

1. A gravitating torus with mass < 0.1 Afc is stable on 
time-scales comparable to the lifetimes of the astrophysical 
objects. 

2. In the general case, the equilibrium cross-section 
of the torus has an oval shape. In equilibrium state, the 
directions of particle movements in the torus cross-section 
are randomized. 

3. The density distribution of particles in the torus 
cross-section obeys Gauss' law with the maximum near to 
the centre of the cross-section. 

In these simulations we assumed that clouds do not 
change their shape for close approaches and collisions. Real 
cloud could significantly change their shape as a result of 
tidal shearing and collisions. Therefore the question arises: 
how can orbiting dusty clouds survive for that long time? 
Krolik & Begelman (1988) considered this problem in detail. 
They argued that tidal shearing can shred the clouds into 
smaller pieces. On the other hand, during collisions the 
material of clouds can be gathered into one giant clump. 



5.1 Obscuration of the central engine by the torus 

Let us find the restrictions on size and number of clouds from 
a requirement of total obscuration in the equatorial plane. 
Assuming an optical depth of clouds in the ultraviolet and 
visible light much greater than unity, we shall consider them 
as an opaque slab with an area of cross-section Tve^R^, where 
the effective radius of a cloud is Rd = eR . The volume per 
one cloud in the case of a circular torus with a homogeneous 
distribution of density has the form 



Vn 



N 



27r^i?Vg 
N ■ 



(37) 



We determine the average number of clouds Ns along the line 
of sight in the equatorial plane of the torus as a product of 
the number density 1/Vn and the volume of a cylinder with 
the cross-sectional area ne^B? and cylinder length equals to 
the diameter of the torus cross-section 



N 9 A'e 
Ns = -y—ne^R^2Rro = — 

Vtorus TTro 



(38) 



However, the cross-section of a torus affected by self- 
gravitation takes an oval shape(Fig.lO) and the distribution 
density of the clouds in its cross-section is inhomogeneous 
(Section 4). Therefore, the expression for the average 
number of clouds along the line of sight will differ from 
(|38|) . In this case, the value of the normalization factor no 
in (|35p can be determined from the following condition: 
the total number of clouds in the torus remains constant 
when the cross-section of the torus is changed owing to self- 
gravitation; that is. 



N = 2Tmo // rfS'exp[-/(r?,C;n))]. 



(39) 



Here the integral is taken over an area of torus cross-section 
S. To obtain the dependence of the number of clouds on 
the angle 9 between the equatorial plane and the line of 
sight we transform into the polar system of coordinates with 
the origin at the central mass: p = rj -\- 1 = rcos9, ( = 
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Figure 10. An average density distribution of particles in an equilibrium cross-section of torus consisting of N=8 192 particles for various 
values of the torus mass and the initial geometrical parameter: a)ro = 0.3, Mtorus = 0.02Mc, b)ro = 0.5, Mtorus = 0.056Mc c) tq = 0.6, 
orus — O.OSA^c, d) vq — 0.7, A'ltorus — 0.109iV/c. The density distributions are plotted in cylindrical coordinates (/y; and the central 
mass is at the point with coordinates (-1, 0). 




Figure 11. Dependence of the number of clouds along the line of 
sight on the equatorial angle 6 calculated by (|40j. The different 
curves correspond to tori with equilibrium cross-sections obtained 
for different values of the torus mass and the initial geometrical 
parameter: l)ro = 0.5, Mtorus = 0.056A{fc 2) rg = 0.6, Mtorus = 
O.OSMo, 3) ro = 0.7, Mtorus = O.IOQM^. 

The total number of clouds in the torus A'^ = 10^ and e = 0.01. 

rain 9. Then the number of clouds along the line of sight as 
a function of angle is 



Ns{e) 



s 



(40) 



Thus, the average number of clouds along the line of sight 
in the equatorial plane of the torus, taking into account the 
inhomogeneous distribution of particles (|40[l . is given by: 



Ns=NsiO) 



2 JJe-f(^-OdS' 

s 



(41) 



For the limiting case of a homogeneous circular torus, the 
integral in the numerator is equal to 2ro and the integral in 
the denominator is equal to the area of a circle, and thus 
the expression (|41[) reduces to (|38|) . The number of clouds 
along the line of sight is an important parameter because 
it characterizes the optical thickness of the torus. Indeed, 
by virtue of the random distribution of clouds in the torus, 
the fraction of non-obscured area Ip, which corresponds to 
the fraction of transmitted radiation of the central engine, 
is determined by a Poisson distribution Ip = prob(k) — 



A^^e~^=/fc! with fc = 0, where k is the frequency of events. 
Consequently, Ip = e^^°. On the other hand, because 
Ip = e~'^^ , the optical depth of the torus in visible hght 
equals the average number of clouds: ry — Ns- Fig. 11 shows 
the number of clouds along the line of sight as a function 
of the equatorial angle 9 for three cases corresponding the 
equilibrium cross-sections in Fig. 10 b,c,d. In order to fit 
the observational data, the number of clouds along the line 
of sight in the equatorial plane of the torus Ns should be 
greater than 5 (Nenkova et al., part 2, 2008). If e = 0.01 
then this condition holds for a torus with total number of 
clouds A'' = 10® and ro > 0.5. The total number of clouds 
decreases with the increasing relative size of clouds because 
N oc e-^. 



5.2 Dynamics of clouds in the torus 

Clouds in the torus are orbiting in the gravitational field of 
the central mass. Owing the fact that the orbits of clouds 
have a scatter in eccentricities and inclinations, they form 
a toroidal structure. Although the orbits are perturbed by 
mutual attraction of the clouds, the average velocity near the 
torus plane of symmetry can be determined as the velocity 
in a circular orbit: 



V ~ 210fcm/s 



Me 



[ipcj 



(42) 



Then the orbital period of a cloud can be estimated as 



-1/2 / ^ N 3/2 



(43) 



For the central mass Mc = IO^Mq the velocity is V = 
150km /s at distance r = 2pc, and the corresponding orbital 
period P ~ SOOOOyear. The inner edge of the torus (rmin) 
is formed by clouds that are moving in elliptical orbits and 
pass through the pericentre, while the outer edge (rmax) 
is formed by clouds that pass through the apocentre. So, 
the cloud velocities in these two regions can be estimated 
approximately as: 



14, 



Vo 



1 + 



ro 



Vtnin — Vo 



- ro 



l + ro' 



(44) 



where Vo = \/ GMc/a and a — {r,nax + rmin)/^ is the semi- 
major axis. For example if rmin = 0.4pc and r^ax = 4pc 
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the maximum velocity of clouds l^max — 330km/ s and the 
minimum velocity Vmin — 60km /s for vq — 0.7. 



6 CONCLUSION 

In this paper the gravitational properties of a torus have 
been investigated. It is shown that the presence of a central 
mass is a necessary condition for the stability of a self- 
gravitating torus. If the mass of the torus is much lower 
than the central mass, such a system can be considered in 
the framework of the problem of test particles motion in 
the central gravitational field. The toroidal distribution of 
particles is achieved as a result of the significant spread of 
inclinations and eccentricities of their orbits. In this case 
the particles move on Keplerian orbits, and we call a torus 
formed in this way as 'Keplerian torus'. It is shown that 
Keplerian torus has a limit on the geometric parameter. 

To investigate the properties of a self-gravitating torus 
we considered the A'^-body problem for a torus located in 
gravitational field of a central mass. As initial conditions 
we used the Keplerian torus. It is shown that the self- 
gravitating torus is stable on time-scale comparable to 
lifetimes of astrophysical objects. The equilibrium cross- 
section of the torus has an oval shape with Gaussian density 
distribution. Until now, the cross-section shape and the 
distribution of particles in a self-gravitating torus were 
not known and therefore were taken as free parameters in 
problems such as radiation transfer in an obscuring torus in 
AGN. We found the dependence of the obscuring efficiency 
as a function of the angle between the line of sight and the 
torus plane of symmetry. 

An account of the dissipation processes is a separate 
issue that is beyond the scope of this paper. The clouds can 
be heated as a result of collisions, thus creating additional 
radiation pressure inside the torus. It is likely that these 
processes can lead to a 'repulsion' of clouds and hence 
additionally influence at the shape of the torus cross-section. 
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